home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 2003 August / MW 8 2003 CD1.iso / Inside Macworld / Product News / gimp-1.2.4.sit / gimp-1.2.4 / plug-ins / ifscompose / ifscompose_utils.c < prev    next >
Encoding:
C/C++ Source or Header  |  2003-02-13  |  21.9 KB  |  925 lines
  1. /* The GIMP -- an image manipulation program
  2.  * Copyright (C) 1995 Spencer Kimball and Peter Mattis
  3.  *
  4.  * IfsCompose is a interface for creating IFS fractals by
  5.  * direct manipulation.
  6.  * Copyright (C) 1997 Owen Taylor
  7.  *
  8.  * it under the terms of the GNU General Public License as published by
  9.  * the Free Software Foundation; either version 2 of the License, or
  10.  * (at your option) any later version.
  11.  *
  12.  * This program is distributed in the hope that it will be useful,
  13.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15.  * GNU General Public License for more details.
  16.  *
  17.  * You should have received a copy of the GNU General Public License
  18.  * along with this program; if not, write to the Free Software
  19.  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  20.  */
  21.  
  22. #include "config.h"
  23.  
  24. #include <stdlib.h>
  25. #include <stdio.h>
  26. #include <string.h>
  27.  
  28. #include <gdk/gdk.h>
  29.  
  30. #include <libgimp/gimp.h>
  31.  
  32. #include "ifscompose.h"
  33.  
  34. typedef struct {
  35.   GdkPoint point;
  36.   gdouble angle;
  37. } SortPoint;
  38.  
  39. /* local functions */
  40. static void 
  41. aff_element_compute_click_boundary(AffElement *elem, int num_elements,
  42.                    gdouble *points_x, gdouble *points_y);
  43. static guchar *
  44. create_brush(IfsComposeVals *ifsvals, gint *brush_size, gdouble *brush_offset);
  45.  
  46. void
  47. aff2_translate(Aff2 *naff, gdouble x, gdouble y)
  48. {
  49.   naff->a11 = 1.0;
  50.   naff->a12 = 0;
  51.   naff->a21 = 0;
  52.   naff->a22 = 1.0;
  53.   naff->b1 = x;
  54.   naff->b2 = y;
  55. }
  56.  
  57. void
  58. aff2_rotate(Aff2 *naff, gdouble theta)
  59. {
  60.   naff->a11 = cos(theta);
  61.   naff->a12 = sin(theta);
  62.   naff->a21 = -naff->a12;
  63.   naff->a22 = naff->a11;
  64.   naff->b1 = 0;
  65.   naff->b2 = 0;
  66. }
  67.  
  68. void
  69. aff2_scale(Aff2 *naff, gdouble s, gint flip)
  70. {
  71.   if (flip)
  72.     naff->a11 = -s;
  73.   else
  74.     naff->a11 = s;
  75.   naff->a12 = 0;
  76.   naff->a21 = 0;
  77.   naff->a22 = s;
  78.   naff->b1 = 0;
  79.   naff->b2 = 0;
  80. }
  81.  
  82. /* Create a unitary transform with given x-y asymmetry and shear */
  83. void 
  84. aff2_distort(Aff2 *naff, gdouble asym, gdouble shear)
  85. {
  86.   naff->a11 = asym;
  87.   naff->a22 = 1/asym;
  88.   naff->a12 = shear;
  89.   naff->a21 = 0;
  90.   naff->b1 = 0;
  91.   naff->b2 = 0;
  92. }
  93.  
  94. /* Find a pure stretch in some directon that brings xo,yo to xn,yn */
  95. void
  96. aff2_compute_stretch(Aff2 *naff,
  97.              gdouble xo, gdouble yo,
  98.              gdouble xn, gdouble yn)
  99. {
  100.   gdouble denom = xo*xn + yo*yn;
  101.   
  102.   if (denom == 0.0)    /* singular */
  103.     {
  104.       naff->a11 = 1.0;
  105.       naff->a12 = 0.0;
  106.       naff->a21 = 0.0;
  107.       naff->a22 = 1.0;
  108.     }
  109.   else
  110.     {
  111.       naff->a11 = (SQR(xn) + SQR(yo))/denom;
  112.       naff->a22 = (SQR(xo) + SQR(yn))/denom;
  113.       naff->a12 = naff->a21 = (xn*yn - xo*yo)/denom;
  114.     }
  115.   
  116.   naff->b1 = 0.0;
  117.   naff->b2 = 0.0;
  118. }
  119.  
  120. void 
  121. aff2_compose(Aff2 *naff, Aff2 *aff1, Aff2 *aff2)
  122. {
  123.   naff->a11 = aff1->a11*aff2->a11 + aff1->a12*aff2->a21;
  124.   naff->a12 = aff1->a11*aff2->a12 + aff1->a12*aff2->a22;
  125.   naff->b1 = aff1->a11*aff2->b1 + aff1->a12*aff2->b2 + aff1->b1;
  126.   naff->a21 = aff1->a21*aff2->a11 + aff1->a22*aff2->a21;
  127.   naff->a22 = aff1->a21*aff2->a12 + aff1->a22*aff2->a22;
  128.   naff->b2 = aff1->a21*aff2->b1 + aff1->a22*aff2->b2 + aff1->b2;
  129. }
  130.  
  131. /* Returns the identity matrix if the original matrix was singular */
  132. void
  133. aff2_invert(Aff2 *naff, Aff2 *aff)
  134. {
  135.   gdouble det = aff->a11*aff->a22 - aff->a12*aff->a21;
  136.   
  137.   if (det==0)
  138.     {
  139.       aff2_scale(naff,1.0,0);
  140.     }
  141.   else
  142.     {
  143.       naff->a11 = aff->a22 / det;
  144.       naff->a22 = aff->a11 / det;
  145.       naff->a21 = - aff->a21 / det;
  146.       naff->a12 = - aff->a12 / det;
  147.       naff->b1 = - naff->a11*aff->b1 - naff->a12*aff->b2;
  148.       naff->b2 = - naff->a21*aff->b1 - naff->a22*aff->b2;
  149.     }
  150. }
  151.  
  152. void 
  153. aff2_apply(Aff2 *aff, gdouble x, gdouble y,
  154.                gdouble *xf, gdouble *yf)
  155. {
  156.   gdouble xt = aff->a11*x + aff->a12*y + aff->b1;
  157.   gdouble yt = aff->a21*x + aff->a22*y + aff->b2;
  158.   *xf = xt;
  159.   *yf = yt;
  160. }
  161.  
  162. /* Find the fixed point of an affine transformation
  163.    (Will return garbage for pure translations) */
  164.  
  165. void
  166. aff2_fixed_point(Aff2 *aff, gdouble *xf, gdouble *yf)
  167. {
  168.   Aff2 t1,t2;
  169.  
  170.   t1.a11 = 1-aff->a11;
  171.   t1.a22 = 1-aff->a22;
  172.   t1.a12 = -aff->a12;
  173.   t1.a21 = -aff->a21;
  174.   t1.b1 = 0;
  175.   t1.b2 = 0;
  176.  
  177.   aff2_invert(&t2,&t1);
  178.   aff2_apply(&t2,aff->b1,aff->b2,xf,yf);
  179. }
  180.  
  181. void 
  182. aff3_apply (Aff3 *t, gdouble x, gdouble y, gdouble z,
  183.         gdouble *xf, gdouble *yf, gdouble *zf)
  184. {
  185.   double xt = t->vals[0][0]*x + t->vals[0][1]*y + t->vals[0][2]*z + t->vals[0][3];
  186.   double yt = t->vals[1][0]*x + t->vals[1][1]*y + t->vals[1][2]*z + t->vals[1][3];
  187.   double zt = t->vals[2][0]*x + t->vals[2][1]*y + t->vals[2][2]*z + t->vals[2][3];
  188.  
  189.   *xf = xt;
  190.   *yf = yt;
  191.   *zf = zt;
  192. }
  193.  
  194. static int
  195. ipolygon_sort_func(const void *a, const void *b)
  196. {
  197.   if (((SortPoint *)a)->angle < ((SortPoint *)b)->angle)
  198.     return -1;
  199.   else if (((SortPoint *)a)->angle > ((SortPoint *)b)->angle)
  200.     return 1;
  201.   else
  202.     return 0;
  203. }
  204.  
  205. /* Return a newly-allocated polygon which is the convex hull
  206.    of the given polygon.
  207.  
  208.    Uses the Graham scan. see
  209.    http://www.cs.curtin.edu.au/units/cg201/notes/node77.html
  210.  
  211.    for a description
  212. */
  213.  
  214. IPolygon *
  215. ipolygon_convex_hull(IPolygon *poly)
  216. {
  217.   gint num_new = poly->npoints;
  218.   GdkPoint *new_points = g_new(GdkPoint,num_new);
  219.   SortPoint *sort_points = g_new(SortPoint,num_new);
  220.   IPolygon *new_poly = g_new(IPolygon,1);
  221.  
  222.   gint i,j;
  223.   gint x1,x2,y1,y2;
  224.   gint lowest;
  225.   GdkPoint lowest_pt;
  226.  
  227.   new_poly->points = new_points;
  228.   if (num_new <= 3)
  229.     {
  230.       memcpy(new_points,poly->points,num_new*sizeof(GdkPoint));
  231.       new_poly->npoints = num_new;
  232.       return new_poly;
  233.     }
  234.  
  235.   /* scan for the lowest point */
  236.   lowest_pt = poly->points[0];
  237.   lowest = 0;
  238.  
  239.   for (i=1;i<num_new;i++)
  240.     if (poly->points[i].y < lowest_pt.y)
  241.       {
  242.     lowest_pt = poly->points[i];
  243.     lowest = i;
  244.       }
  245.  
  246.   /* sort by angle from lowest point */
  247.   
  248.   for (i=0,j=0;i<num_new;i++,j++)
  249.     {
  250.       if (i==lowest)
  251.     j--;
  252.       else
  253.     {
  254.       gdouble dy = poly->points[i].y - lowest_pt.y;
  255.       gdouble dx = poly->points[i].x - lowest_pt.x;
  256.  
  257.       if (dy==0 && dx==0)
  258.         {
  259.           j--;
  260.           num_new--;
  261.           continue;
  262.         }
  263.       sort_points[j].point = poly->points[i];
  264.       sort_points[j].angle = atan2(dy,dx);
  265.     }
  266.     }
  267.  
  268.   qsort(sort_points,num_new-1,sizeof(SortPoint),ipolygon_sort_func);
  269.  
  270.   /* now ensure that all turns as we trace the perimiter are
  271.      counter-clockwise */
  272.  
  273.   new_points[0] = lowest_pt;
  274.   new_points[1] = sort_points[0].point;
  275.   x1 = new_points[1].x - new_points[0].x;
  276.   y1 = new_points[1].y - new_points[0].y;
  277.  
  278.   for (i=1,j=2;j<num_new;i++,j++)
  279.     {
  280.       x2 = sort_points[i].point.x - new_points[j-1].x;
  281.       y2 = sort_points[i].point.y - new_points[j-1].y;
  282.  
  283.       if (x2==0 && y2==0)
  284.     {
  285.       num_new--;
  286.       j--;
  287.       continue;
  288.     }
  289.       
  290.       while (x1*y2 - x2*y1 < 0)    /* clockwise rotation */
  291.     {
  292.       num_new--;
  293.       j--;
  294.       x1 = new_points[j-1].x - new_points[j-2].x;
  295.       y1 = new_points[j-1].y - new_points[j-2].y;
  296.       x2 = sort_points[i].point.x - new_points[j-1].x;
  297.       y2 = sort_points[i].point.y - new_points[j-1].y;
  298.     }
  299.       new_points[j] = sort_points[i].point;
  300.       x1 = x2;
  301.       y1 = y2;
  302.     }
  303.  
  304.   g_free(sort_points);
  305.  
  306.   new_poly->npoints = num_new;
  307.   return new_poly;
  308. }
  309.  
  310. /* Determines whether a specified point is in the given polygon.
  311.    Based on
  312.  
  313.    inpoly.c by Bob Stein and Craig Yap.
  314.  
  315.    (Linux Journal, Issue 35 (March 1997), p 68)
  316.    */
  317.  
  318. gint
  319. ipolygon_contains(IPolygon *poly, gint xt, gint yt)
  320. {
  321.   gint xnew, ynew;
  322.   gint xold, yold;
  323.   gint x1,y1;
  324.   gint x2,y2;
  325.  
  326.   gint i;
  327.   gint inside = 0;
  328.  
  329.   if (poly->npoints < 3)
  330.     return 0;
  331.  
  332.   xold=poly->points[poly->npoints-1].x;
  333.   yold=poly->points[poly->npoints-1].y;
  334.   for (i=0;i<poly->npoints;i++)
  335.     {
  336.       xnew = poly->points[i].x;
  337.       ynew = poly->points[i].y;
  338.       if (xnew > xold) 
  339.     {
  340.       x1 = xold;
  341.       x2 = xnew;
  342.       y1 = yold;
  343.       y2 = ynew;
  344.     }
  345.       else
  346.     {
  347.       x1 = xnew;
  348.       x2 = xold;
  349.       y1 = ynew;
  350.       y2 = yold;
  351.     }
  352.       if ((xnew < xt) == (xt <= xold) &&
  353.       (yt - y1)*(x2 - x1) < (y2 - y1)*(xt - x1))
  354.     inside = !inside;
  355.       xold = xnew;
  356.       yold = ynew;
  357.     }
  358.   return inside;
  359. }
  360.  
  361. void
  362. aff_element_compute_color_trans(AffElement *elem)
  363. {
  364.   int i,j;
  365.  
  366.   if (elem->v.simple_color)
  367.     {
  368.       gdouble mag2 = 0;
  369.       for (i=0;i<3;i++)
  370.     mag2 += SQR(elem->v.target_color.vals[i]);
  371.  
  372.       /* For mag2 == 0, the transformation blows up in general
  373.      but is well defined for hue_scale == value_scale, so
  374.      we assume that special case. */
  375.       if (mag2 == 0)
  376.     for (i=0;i<3;i++)
  377.       {
  378.         for (j=0;j<4;j++)
  379.           elem->color_trans.vals[i][j] = 0.0;
  380.         elem->color_trans.vals[i][i] = elem->v.hue_scale;
  381.       }
  382.       else
  383.     for (i=0;i<3;i++)
  384.       {
  385.         for (j=0;j<3;j++)
  386.           {
  387.         elem->color_trans.vals[i][j] = elem->v.target_color.vals[i]
  388.           / mag2 * (elem->v.value_scale - elem->v.hue_scale);
  389.         if (i==j)
  390.           elem->color_trans.vals[i][j] += elem->v.hue_scale;
  391.           }
  392.         elem->color_trans.vals[i][3] = 
  393.           (1-elem->v.value_scale)*elem->v.target_color.vals[i];
  394.       }
  395.       aff3_apply(&elem->color_trans,1.0,0.0,0.0,&elem->v.red_color.vals[0],
  396.          &elem->v.red_color.vals[1],&elem->v.red_color.vals[2]);
  397.       aff3_apply(&elem->color_trans,0.0,1.0,0.0,&elem->v.green_color.vals[0],
  398.          &elem->v.green_color.vals[1],&elem->v.green_color.vals[2]);
  399.       aff3_apply(&elem->color_trans,0.0,0.0,1.0,&elem->v.blue_color.vals[0],
  400.          &elem->v.blue_color.vals[1],&elem->v.blue_color.vals[2]);
  401.       aff3_apply(&elem->color_trans,0.0,0.0,0.0,&elem->v.black_color.vals[0],
  402.          &elem->v.black_color.vals[1],&elem->v.black_color.vals[2]);
  403.     }
  404.   else
  405.     {
  406.       for (i=0;i<3;i++)
  407.     elem->color_trans.vals[i][0] = elem->v.red_color.vals[i]
  408.       - elem->v.black_color.vals[i];
  409.       for (i=0;i<3;i++)
  410.     elem->color_trans.vals[i][1] = elem->v.green_color.vals[i]
  411.       - elem->v.black_color.vals[i];
  412.       for (i=0;i<3;i++)
  413.     elem->color_trans.vals[i][2] = elem->v.blue_color.vals[i]
  414.       - elem->v.black_color.vals[i];
  415.       for (i=0;i<3;i++)
  416.     elem->color_trans.vals[i][3] = elem->v.black_color.vals[i];
  417.     }
  418. }
  419.  
  420. void
  421. aff_element_compute_trans(AffElement *elem, gdouble width, gdouble height,
  422.               gdouble center_x, gdouble center_y)
  423. {
  424.   Aff2 t1, t2, t3;
  425.  
  426.   /* create the rotation, scaling and shearing part of the transform */
  427.   aff2_distort(&t1, elem->v.asym, elem->v.shear);
  428.   aff2_scale(&t2, elem->v.scale, elem->v.flip);
  429.   aff2_compose(&t3, &t2, &t1);
  430.   aff2_rotate(&t2, elem->v.theta);
  431.   aff2_compose(&t1, &t2, &t3);
  432.  
  433.   /* now create the translational part */
  434.   aff2_translate(&t2, -center_x*width, -center_y*width);
  435.   aff2_compose(&t3, &t1, &t2);
  436.   aff2_translate(&t2, elem->v.x*width, elem->v.y*width);
  437.   aff2_compose(&elem->trans, &t2, &t3);
  438. }
  439.  
  440. void
  441. aff_element_decompose_trans(AffElement *elem, Aff2 *aff, gdouble width,
  442.                 gdouble height, gdouble center_x,
  443.                 gdouble center_y)
  444. {
  445.   Aff2 t1,t2;
  446.   gdouble det,scale,sign;
  447.  
  448.   /* pull of the translational parts */
  449.   aff2_translate(&t1,center_x*width,center_y*width);
  450.   aff2_compose(&t2,aff,&t1);
  451.   
  452.   elem->v.x = t2.b1 / width;
  453.   elem->v.y = t2.b2 / width;
  454.  
  455.   det = t2.a11*t2.a22 - t2.a12*t2.a21;
  456.  
  457.   if (det == 0.0)
  458.     {
  459.       elem->v.scale = 0.0;
  460.       elem->v.theta = 0.0;
  461.       elem->v.asym = 1.0;
  462.       elem->v.shear = 0.0;
  463.       elem->v.flip = 0;
  464.     }
  465.   else
  466.     {
  467.       if (det >= 0)
  468.     {
  469.       scale = elem->v.scale = sqrt(det);
  470.       sign = 1;
  471.       elem->v.flip = 0;
  472.     }
  473.       else
  474.     {
  475.       scale = elem->v.scale = sqrt(-det);
  476.       sign = -1;
  477.       elem->v.flip = 1;
  478.     }
  479.  
  480.       elem->v.theta = atan2(-t2.a21,t2.a11);
  481.  
  482.       if (cos(elem->v.theta) == 0.0)
  483.     {
  484.       elem->v.asym = - t2.a21 / scale / sin(elem->v.theta);
  485.       elem->v.shear = - sign * t2.a22 / scale / sin(elem->v.theta);
  486.     }
  487.       else
  488.     {
  489.       elem->v.asym = sign * t2.a11 / scale / cos(elem->v.theta);
  490.       elem->v.shear = sign * 
  491.         (t2.a12/scale - sin(elem->v.theta)/elem->v.asym)
  492.         / cos(elem->v.theta);
  493.     }
  494.     }
  495. }
  496.  
  497. static void 
  498. aff_element_compute_click_boundary(AffElement *elem, int num_elements,
  499.                    gdouble *points_x, gdouble *points_y)
  500. {
  501.   gint i;
  502.   gdouble xtot = 0;
  503.   gdouble ytot = 0;
  504.   gdouble xc, yc;
  505.   gdouble theta;
  506.   gdouble sth,cth;        /* sin(theta), cos(theta) */
  507.   gdouble axis1,axis2;
  508.   gdouble axis1max, axis2max, axis1min, axis2min;
  509.  
  510.   /* compute the center of mass of the points */
  511.   for (i=0; i<num_elements; i++)
  512.     {
  513.       xtot += points_x[i];
  514.       ytot += points_y[i];
  515.     }
  516.   xc = xtot/num_elements;
  517.   yc = ytot/num_elements;
  518.  
  519.   /* compute the sum of the (x+iy)^2, and take half the the resulting
  520.      angle (xtot+iytot = A*exp(2i*theta)), to get an average direction */
  521.  
  522.   xtot = 0;
  523.   ytot = 0;
  524.   for (i=0; i<num_elements; i++)
  525.     {
  526.       xtot += SQR(points_x[i]-xc)-SQR(points_y[i]-yc);
  527.       ytot += 2*(points_x[i]-xc)*(points_y[i]-yc);
  528.     }
  529.   theta = 0.5*atan2(ytot,xtot);
  530.   sth = sin(theta);
  531.   cth = cos(theta);
  532.  
  533.   /* compute the minimum rectangle at angle theta that bounds the points,
  534.      1/2 side lenghs left in axis1, axis2, center in xc, yc */
  535.  
  536.   axis1max = axis1min = 0.0;
  537.   axis2max = axis2min = 0.0;
  538.   for (i=0; i<num_elements; i++)
  539.     {
  540.       gdouble proj1 = (points_x[i]-xc)*cth + (points_y[i]-yc)*sth;
  541.       gdouble proj2 = -(points_x[i]-xc)*sth + (points_y[i]-yc)*cth;
  542.       if (proj1 < axis1min)
  543.     axis1min = proj1;
  544.       if (proj1 > axis1max)
  545.     axis1max = proj1;
  546.       if (proj2 < axis2min)
  547.     axis2min = proj2;
  548.       if (proj2 > axis2max)
  549.     axis2max = proj2;
  550.     }
  551.   axis1 = 0.5*(axis1max - axis1min);
  552.   axis2 = 0.5*(axis2max - axis2min);
  553.   xc += 0.5*((axis1max + axis1min)*cth - (axis2max+axis2min)*sth);
  554.   yc += 0.5*((axis1max + axis1min)*sth + (axis2max+axis2min)*cth);
  555.  
  556.   /* if the the rectangle is less than 10 pixels in any dimension,
  557.      make it click_boundary, otherwise set click_boundary = draw_boundary */
  558.  
  559.   if (axis1 < 8.0 || axis2 < 8.0)
  560.     {
  561.       GdkPoint *points = g_new(GdkPoint,4);
  562.  
  563.       elem->click_boundary = g_new(IPolygon,1);
  564.       elem->click_boundary->points = points;
  565.       elem->click_boundary->npoints = 4;
  566.  
  567.       if (axis1 < 8.0) axis1 = 8.0;
  568.       if (axis2 < 8.0) axis2 = 8.0;
  569.       
  570.       points[0].x = xc + axis1*cth - axis2*sth;
  571.       points[0].y = yc + axis1*sth + axis2*cth;
  572.       points[1].x = xc - axis1*cth - axis2*sth;
  573.       points[1].y = yc - axis1*sth + axis2*cth;
  574.       points[2].x = xc - axis1*cth + axis2*sth;
  575.       points[2].y = yc - axis1*sth - axis2*cth;
  576.       points[3].x = xc + axis1*cth + axis2*sth;
  577.       points[3].y = yc + axis1*sth - axis2*cth;
  578.     }
  579.   else
  580.     elem->click_boundary = elem->draw_boundary;
  581. }
  582.  
  583. void 
  584. aff_element_compute_boundary(AffElement *elem, gint width,
  585.                  gint height,
  586.                  AffElement **elements, 
  587.                  int num_elements)
  588. {
  589.   int i;
  590.   IPolygon tmp_poly;
  591.   gdouble *points_x;
  592.   gdouble *points_y;
  593.  
  594.   if (elem->click_boundary && elem->click_boundary != elem->draw_boundary)
  595.     g_free(elem->click_boundary);
  596.   if (elem->draw_boundary)
  597.     g_free(elem->draw_boundary);
  598.  
  599.   tmp_poly.npoints = num_elements;
  600.   tmp_poly.points = g_new(GdkPoint,num_elements);
  601.   points_x = g_new(gdouble,num_elements);
  602.   points_y = g_new(gdouble,num_elements);
  603.   
  604.   for (i=0;i<num_elements;i++)
  605.     {
  606.       aff2_apply(&elem->trans,elements[i]->v.x*width,elements[i]->v.y*width,
  607.          &points_x[i],&points_y[i]);
  608.       tmp_poly.points[i].x = (gint)points_x[i];
  609.       tmp_poly.points[i].y = (gint)points_y[i];
  610.     }
  611.  
  612.   elem->draw_boundary = ipolygon_convex_hull(&tmp_poly);
  613.   aff_element_compute_click_boundary(elem,num_elements,points_x,points_y);
  614.  
  615.   g_free(tmp_poly.points);
  616. }
  617.  
  618. void 
  619. aff_element_draw(AffElement *elem, gint selected,
  620.          gint width, gint height, 
  621.          GdkDrawable *win,
  622.          GdkGC *normal_gc,GdkGC *selected_gc,
  623.          GdkFont *font)
  624. {
  625.   GdkGC *gc;
  626.   gint string_width = gdk_string_width (font,elem->name);
  627.   gint string_height = font->ascent  + font->descent + 2;
  628.  
  629.   if (selected)
  630.     gc = selected_gc;
  631.   else
  632.     gc = normal_gc;
  633.  
  634.   gdk_draw_string(win,font,gc,
  635.           elem->v.x*width-string_width/2,
  636.           elem->v.y*width+string_height/2,elem->name);
  637.  
  638.   if (elem->click_boundary != elem->draw_boundary)
  639.     gdk_draw_polygon(win,normal_gc,FALSE,elem->click_boundary->points,
  640.              elem->click_boundary->npoints);
  641.  
  642.   gdk_draw_polygon(win,gc,FALSE,elem->draw_boundary->points,
  643.            elem->draw_boundary->npoints);
  644. }
  645.  
  646. AffElement *
  647. aff_element_new(gdouble x, gdouble y, IfsColor color, gint count)
  648. {
  649.   AffElement *elem = g_new(AffElement, 1);
  650.   char buffer[16];
  651.  
  652.   elem->v.x = x;
  653.   elem->v.y = y;
  654.   elem->v.theta = 0.0;
  655.   elem->v.scale = 0.5;
  656.   elem->v.asym = 1.0;
  657.   elem->v.shear = 0.0;
  658.   elem->v.flip = 0;
  659.   
  660.   elem->v.red_color = color;
  661.   elem->v.blue_color = color;
  662.   elem->v.green_color = color;
  663.   elem->v.black_color = color;
  664.  
  665.   elem->v.target_color = color;
  666.   elem->v.hue_scale = 0.5;
  667.   elem->v.value_scale = 0.5;
  668.  
  669.   elem->v.simple_color = TRUE;
  670.  
  671.   elem->draw_boundary = NULL;
  672.   elem->click_boundary = NULL;
  673.  
  674.   aff_element_compute_color_trans(elem);
  675.  
  676.   elem->v.prob = 1.0;
  677.  
  678.   sprintf(buffer,"%d",count);
  679.   elem->name = g_strdup(buffer);
  680.   
  681.   return elem;
  682. }
  683.  
  684. void
  685. aff_element_free(AffElement *elem)
  686. {
  687.   if (elem->click_boundary != elem->draw_boundary)
  688.     g_free(elem->click_boundary);
  689.   g_free(elem->draw_boundary);
  690.   g_free(elem);
  691. }
  692.  
  693. #ifdef DEBUG_BRUSH
  694. static brush_chars[] = {' ',':','*','@'};
  695. #endif
  696.  
  697. static guchar *
  698. create_brush(IfsComposeVals *ifsvals, gint *brush_size, gdouble *brush_offset)
  699. {
  700.   gint i,j;
  701.   gint ii,jj;
  702.   guchar *brush;
  703. #ifdef DEBUG_BRUSH
  704.   gdouble totpix = 0.0;
  705. #endif
  706.   
  707.   gdouble radius = ifsvals->radius * ifsvals->subdivide;
  708.   *brush_size = ceil(2*radius);
  709.   *brush_offset = 0.5 * (*brush_size-1);
  710.  
  711.   brush = g_new(guchar,SQR(*brush_size));
  712.   
  713.   for (i=0;i<*brush_size;i++)
  714.     {
  715.       for (j=0;j<*brush_size;j++)
  716.     {
  717.       gdouble d = sqrt(SQR(i-*brush_offset)+SQR(j-*brush_offset));
  718.       gdouble pixel = 0.0;
  719.  
  720.       if (d-0.5*sqrt(2) > radius)
  721.         pixel = 0.0;
  722.       else if (d+0.5*sqrt(2) < radius)
  723.         pixel = 1.0;
  724.       else
  725.         for (ii=0;ii<10;ii++)
  726.           for (jj=0;jj<10;jj++)
  727.         {
  728.           d = sqrt(SQR(i-*brush_offset+ii*0.1-0.45)
  729.                +SQR(j-*brush_offset+jj*0.1-0.45));
  730.           pixel += (d<radius)/100.0;
  731.         }
  732.       brush[i**brush_size+j] = 255.999*pixel;
  733. #ifdef DEBUG_BRUSH
  734.       putchar(brush_chars[(int)(pixel*3.999)]);
  735.       totpix += pixel;
  736. #endif /* DEBUG_BRUSH */
  737.     }
  738. #ifdef DEBUG_BRUSH
  739.       putchar('\n');
  740. #endif /* DEBUG_BRUSH */
  741.     }
  742. #ifdef DEBUG_BRUSH
  743.   printf("Brush total / area = %f\n",totpix/SQR(ifsvals->subdivide));
  744. #endif /* DEBUG_BRUSH */
  745.   return brush;
  746. }
  747.  
  748. void
  749. ifs_render(AffElement **elements, gint num_elements,
  750.        gint width, gint height, gint nsteps,
  751.        IfsComposeVals *vals, gint band_y, gint band_height,
  752.        guchar *data, guchar *mask, guchar *nhits, gint preview)
  753. {
  754.   gint i,k;
  755.   gdouble x,y;
  756.   gdouble r,g,b;
  757.   gint ri,gi,bi;
  758.   gint p0,psum;
  759.   gdouble pt;
  760.   guchar *ptr;
  761.   gint *prob;
  762.   gdouble *fprob;
  763.   gint subdivide;
  764.  
  765.   guchar *brush = NULL;
  766.   gint brush_size;
  767.   gdouble brush_offset;
  768.  
  769.   if (preview)
  770.     subdivide = 1;
  771.   else
  772.     subdivide = vals->subdivide;
  773.   
  774.   /* compute the probabilities and transforms */
  775.   fprob = g_new(gdouble,num_elements);
  776.   prob = g_new(gint,num_elements);
  777.   pt = 0.0;
  778.   for (i=0;i<num_elements;i++)
  779.     {
  780.       aff_element_compute_trans(elements[i],width*subdivide,height*subdivide,
  781.                 vals->center_x, vals->center_y);
  782.       fprob[i] = fabs(
  783.     elements[i]->trans.a11 * elements[i]->trans.a22
  784.     - elements[i]->trans.a12 * elements[i]->trans.a21);
  785.       /* As a heuristic, if the determinant is really small, it's
  786.      probably a line element, so increase the probability so
  787.      it gets rendered */
  788.       /* FIXME: figure out what 0.01 really should be */
  789.       if (fprob[i] < 0.01) fprob[i] = 0.01;
  790.       fprob[i] *= elements[i]->v.prob;
  791.       
  792.       pt += fprob[i];
  793.     }
  794.  
  795.   psum = 0;
  796.   for (i=0;i<num_elements;i++)
  797.     {
  798.       psum += G_MAXRAND * (fprob[i]/pt);
  799.       prob[i] = psum;
  800.     }
  801.   prob[i-1] = G_MAXRAND;    /* make sure we don't get bitten
  802.                    by roundoff*/
  803.   /* create the brush */
  804.   if (!preview)
  805.     brush = create_brush(vals,&brush_size,&brush_offset);
  806.  
  807.   x = y = 0;
  808.   r = g = b = 0;
  809.  
  810.   /* now run the iteration */
  811.   for (i=0;i<nsteps;i++) 
  812.     {
  813.       if (!preview && !(i % 5000))
  814.     gimp_progress_update ((gdouble) i / (gdouble) nsteps);
  815. #ifdef G_OS_WIN32
  816.       p0 = RAND_FUNC();
  817. #else
  818.       p0 = rand();
  819. #endif
  820.       k=0;
  821.       while (p0 > prob[k])
  822.     k++;
  823.  
  824.       aff2_apply(&elements[k]->trans,x,y,&x,&y);
  825.       aff3_apply(&elements[k]->color_trans,r,g,b,&r,&g,&b);
  826.       if (i<50) continue;
  827.  
  828.       ri= (gint)(255.999*r);
  829.       gi = (gint)(255.999*g);
  830.       bi = (gint)(255.999*b);
  831.       
  832.       if (preview)
  833.     {
  834.       if ((x<width) && (y<(band_y+band_height)) && 
  835.           (x >= 0) && (y >= band_y) &&
  836.           (ri >= 0) && (ri < 256) &&
  837.           (gi >= 0) && (gi < 256) &&
  838.           (bi >= 0) && (bi < 256))
  839.         {
  840.           ptr = data + 3 * (((gint)(y-band_y))*width + (gint)x);
  841.           *ptr++ = ri;
  842.           *ptr++ = gi;
  843.           *ptr = bi;
  844.         }
  845.     }
  846.       else
  847.     if ((ri >= 0) && (ri < 256) &&
  848.         (gi >= 0) && (gi < 256) &&
  849.         (bi >= 0) && (bi < 256))
  850.       {
  851.         guint m_old;
  852.         guint m_new;
  853.         guint m_pix;
  854.         guint n_hits;
  855.         guint old_scale;
  856.         guint pix_scale;
  857.         
  858.         gint index;
  859.         gint ii,jj;
  860.         gint jj0 = floor(y-brush_offset-band_y*subdivide);
  861.         gint ii0 = floor(x-brush_offset);
  862.         gint jjmax,iimax;
  863.         gint jjmin = 0;
  864.         gint iimin = 0;
  865.  
  866.         if (ii0 < 0)
  867.           iimin = - ii0;
  868.         else
  869.           iimin = 0;
  870.         
  871.         if (jj0 < 0)
  872.           jjmin = - jj0;
  873.         else
  874.           jjmin = 0;
  875.         
  876.         if (jj0+brush_size >= subdivide*band_height) 
  877.           jjmax = subdivide*band_height - jj0;
  878.         else
  879.           jjmax = brush_size;
  880.  
  881.         if (ii0+brush_size >= subdivide*width) 
  882.           iimax = subdivide*width - ii0;
  883.         else
  884.           iimax = brush_size;
  885.  
  886.         for (jj=jjmin;jj<jjmax;jj++)
  887.           for (ii=iimin;ii<iimax;ii++)
  888.         {
  889.           index = (jj0+jj)*width*subdivide + ii0 + ii;
  890.           n_hits = nhits[index];
  891.           if (n_hits == 255)
  892.             continue;
  893.           
  894.           m_pix = brush[jj*brush_size+ii];
  895.           if (!m_pix)
  896.             continue;
  897.           nhits[index] = ++n_hits;
  898.           m_old = mask[index];
  899.           m_new = m_old + m_pix - m_old*m_pix/255;
  900.           mask[index] = m_new;
  901.  
  902.           /* relative probability that old colored pixel is on top */
  903.           old_scale = m_old*(255*n_hits-m_pix);
  904.           /* relative probability that new colored pixel is on top */
  905.           pix_scale = m_pix*((255-m_old)*n_hits+m_old);
  906.             
  907.           ptr = data + 3*index;
  908.           *ptr = ( old_scale * (*ptr) + pix_scale * ri ) /
  909.             ( old_scale + pix_scale );
  910.           ptr++;
  911.           *ptr = ( old_scale * (*ptr) + pix_scale * gi ) /
  912.             ( old_scale + pix_scale );
  913.           ptr++;
  914.           *ptr = ( old_scale * (*ptr) + pix_scale * bi ) /
  915.             ( old_scale + pix_scale );
  916.         }
  917.       }
  918.     } /* main iteration */
  919.  
  920.   if (brush)
  921.     g_free(brush);
  922.   g_free(prob);
  923.   g_free(fprob);
  924. }
  925.